function [f] = f_density_power(a, q_bar, qmin, betaq, sig, q_up)
%LOGLIKELIHOOD
% inputs:  b,as(asset),betaqT,c,sig,cs,c0,c1,tauT,kT
% outputs: logLikelihood
%   see the model file

c = (betaq-1)/qmin^(1-betaq);

if false
    a = 16;
    b = 1;
    betaq = 3
    c = 0.7
    sig = 0.02
    tau = 0.01
    k = 0.0004
    beta = 87
    a = apre
    b = apost
    bpre
    bpost
    
end

Phi  = @(z) normcdf(z);
phi  = @(z) normpdf(z);


% Making sure to have column vectors
a = a(:);
gamma = sig * (betaq - 1);
delta = 1-betaq;

if q_up ~= q_bar

    my1 = c .* a.^(-betaq) .* exp(0.5*delta^2*sig^2);
    PHI_down = Phi(sig^-1 .*(log(a)-log(q_bar))-gamma);
    PHI_up   = Phi(sig^-1 *(log(a)-log(q_up))-gamma);
    phi_down =  phi(sig.^-1 .* (log(a)-log(q_bar)));
    my2  = (sig .*a).^-1 * (c/delta) .* (q_up^delta - q_bar^delta);
   
    f = my1.*(1-PHI_down+PHI_up) + my2.*phi_down;
else 
  f = c .* a.^(-betaq) .* exp(0.5.*delta^2.*sig^2) ; %.* Phi(sig.^-1 .* (log(a)-log(qmin)) - (betaq-1).*sig)
end

